1. Field of the Invention
The present invention relates to device simulation for simulating electrical characteristics of a semiconductor device. More specifically, the present invention relates to reduction in calculation time and improvement in accuracy of device simulation for predicting impact ionization and electrical characteristics incidental thereto.
2. Description of the Background Art
Generally, in manufacturing semiconductor devices, device simulation by a computer is widely performed prior to actual trial manufacture. Device simulation is for presuming behavior of carriers in the device and for calculating device characteristics, based on the physical geometry and impurity distribution of a given device.
As the semiconductor device is increasingly miniaturized, electric field in the device increases, and therefore deterioration of the device caused by carriers (hot carriers) accelerated by the high electric field comes to be a serious problem. Since information related to behavior of carriers in the device can be obtained by device simulation, such simulation is an effective method of obtaining knowledge related to device deterioration. Impact ionization refers to a phenomenon in which electron-hole pairs are generated from crystal lattice because of high energy carriers, which is a representative of phenomena caused by hot carriers. As impact ionization is deeply related to deterioration of semiconductor devices, impact ionization must be accurately represented in order to correctly perform simulation of miniaturized semiconductor devices.
Device simulation includes the following methods of analysis (A1) to (A5).
(A1) Drift-diffusion model PA1 (A2) Hydrodynamic model or energy transport model PA1 (A3) Method of analysis according to distribution function expansion PA1 (A4) Method of analysis according to distribution function repetition PA1 (A5) Monte Carlo method
Of these methods of analysis (A1) to (A5), latter ones have higher accuracies. Methods of analysis (A3), (A4) and (A5) include solution of Boltzmann transport equation, and therefore distribution function of the carrier can be calculated.
In the method of analysis (A1), i.e., drift-diffusion model, potential and carrier concentration can be calculated by combining the following Poisson's equation (11) and current continuity equation (12). EQU .gradient..multidot.D=.rho. (11) EQU .gradient..multidot.J=q.multidot.(.differential.n/.differential.t).sub.GR( 12)
In equation (11), the reference character D represents electric flux density, and .rho. represents density of electric charge. The electric flux density D can be represented by the following equation (11a). EQU D=-.epsilon..gradient..psi. (11a)
where .psi. represents potential and .epsilon. represents dielectric constant. The density of electric charge .rho. in equation (11) is determined from hole density, electron density, doner density and acceptor density.
In equation (12), the reference character J represents current density, q denotes a unit charge, n denotes carrier density and (.differential.n/.differential.t).sub.GR denotes net generation rate of carriers. The current density and net generation rate of carriers can be represented by the following equations (12a) and (12b), respectively. EQU J=-n q.mu..gradient..psi.+q D.sub.d .gradient.n (12a) EQU (.differential.n/.omega.t).sub.GR =.alpha..vertline.J/q.vertline.+(.differential.n/.differential.t) .sub.GRO( 12b)
where .mu. represents mobility of the carriers, D.sub.d represents diffusion coefficient, and (.differential.n/.differential.t).sub.GRO denotes net generation rate of carriers caused by mechanisms other than impact ionization.
By using the method of analysis (A1), that is, drift-diffusion model, solution can be obtained in a shortest period of calculation time as compared with other methods of analysis (A2) to (A5). However, by the method of analysis (A1), information related to the energy of carriers cannot be obtained, and therefore it is difficult to accurately predict behavior of hot carriers.
In the method of analysis (A2), that is, hydrodynamic model or energy transport model, potential, carrier density and mean energy of carriers (or carrier temperature) can be calculated by combining Poisson's equation (11), current continuity equation (12) and equation of conservation of energy (13) below. EQU .gradient..multidot.S=E.multidot.J+n(.differential.w/.differential.t).sub.c oll ( 13)
where S denotes energy flow density, E denotes electric field, w denotes energy and (.differential.w/.differential.t).sub.coll denotes change rate of energy caused by scattering of carriers. In the hydrodynamic model or energy transport model of the method of analysis (A2), mean energy of carriers can be obtained, and therefore as compared with the drift-diffusion model of method of analysis (A1), accuracy in predicting hot carrier phenomenon is improved. However, in the method of analysis (A2), the equation of conservation of energy which is highly non-linear must be solved, so that calculation time is increased as compared with the model of the method (A1), and it takes longer until convergence. In the models of methods (A1) and (A2) of analysis, only the mean values (carrier concentration, carrier mean energy) with respect to the carriers are addressed.
In the method of analysis (A3) according to distribution function expansion, carrier distribution in momentum space is represented as expanded in an appropriate function, and coefficients in the function of expansion are calculated for various point of the device to be analyzed. In this method according to distribution function expansion, carrier distribution in the momentum space can be known. However, in order to know accurate carrier distribution, function to be used for expansion and the number of terms expanded must be carefully selected.
In the method of analysis (A4) according to distribution function repetition, momentum space is made discrete, and carrier distribution is calculated by repetitively changing carrier distribution at each discrete point by acceleration by electric field and scattering according to various mechanisms of scattering.
In the method of analysis (A5), that is, Monte Carlo method, individual carrier is moved in the momentum space and in real space by acceleration by the electric field and scattering according to various mechanisms of scattering, and the results of these movement are statistically processed to know the behavior of carriers in the device. Since movement of individual carrier can be known in the Monte Carlo method, hot carrier phenomenon can be analyzed with highest accuracy. Meanwhile, in order to obtain a statistically meaningful solution, it is necessary to calculate movements of large number of carriers for a long period of time, and therefore the necessary calculation time is longer in the order of two digits, than the time necessary in the methods (A1) and (A2) of analysis.
Generally speaking, methods of analysis (A3), (A4) and (A5) provide more accurate solution than the methods (A1) and (A2). However the former methods require longer period of time for calculation.
As already mentioned, impact ionization is caused by hot carriers. Carriers generated by this phenomenon may cause substrate current or gate current. Since carriers injected to the gate is closely related to device deterioration, accurate estimation of impact ionization is very important in view of reliability over long period of time of the device.
Impact ionization is caused by carriers having high energy. Therefore, it is necessary to know energy distribution of carriers in order to calculate possibility of impact ionization. By calculation according to distribution function expansion model or Monte Carlo model, it is possible to know energy distribution of carriers, and therefore impact ionization coefficient (number of impact ionization caused by one carrier while the carrier travels a unit length) can be accurately known. Meanwhile, in the drift-diffusion model, hydrodynamic model or energy transport model, energy distribution of carriers is not known, and therefore some approximation is necessary for calculating impact ionization coefficient.
Device simulation widely performed at present is based on the drift-diffusion model of the method of analysis (A1). However, according to this model, it is impossible to know the energy distribution of carriers. The method of calculating impact ionization coefficient .alpha. using drift-diffusion model includes a method in which .alpha. is represented as a function of magnitude of an electric field of the position to be calculated, and a method of calculating effective electric field by tracing back the current or electric flow line in a direction reverse to the direction of travel of the carriers. The former method provides accurate solution only when the electric field is uniform. Since electric field varies spatially in actual semiconductor devices, accurate solution cannot be obtained by this method. The latter method requires long period of time in calculating impact ionization coefficient, since current or electric flow line must be traced back every time the impact ionization coefficient is to be calculated.
FIG. 5 shows comparison between impact ionization coefficient calculated as a function solely of the magnitude of the electric field, and the impact ionization coefficient calculated in accordance with Monte Carlo method, which is the model providing highest accuracy. In the graph of FIG. 5, the abscissa represents an inverse (cm/MV) of the electric field E, and the ordinate represents impact ionization coefficient .alpha. (cm-1). In FIG. 5, it is assumed that the electric field E is represented by the following exponential function. EQU E=E.sub.0 .multidot.exp(x/.lambda.) (14)
where .lambda. represents characteristic length. Under the electric field thus assumed, lines 5a, 5b, 5d and 5e in the graph represent impact ionization coefficients calculated in accordance with Monte Carlo method when .lambda. is -0.1 .mu.m, -0.2 .mu.m. 0.2 .mu.m and 0.1 .mu.m, respectively. Meanwhile, line 5C represents impact ionization coefficient calculated in a spatially uniform electric field.
It can be seen from FIG. 5 that when the value .lambda. is positive, that is, when electric field is increasing, the impact ionization coefficient (line 5C) calculated based only on the magnitude of the electric field is estimated to be larger than the correct value (as represented by lines 5D and 5E) calculated in accordance with the Monte Carlo method. Further, it can be seen that there is different proportional relation held between 1n.alpha. and 1/E, dependent on the magnitude of the characteristic length .lambda.. Therefore, it is understood that calculation of impact ionization coefficient based only on the magnitude of electric field is not very accurate.
Therefore, in order to calculate impact ionization coefficient accurately, it is necessary to use a method of tracing current or electric flow line, the method according to energy transport model, the method of solving Boltzmann transport equation or the like in order to take into account the variation in the electric field. However, as compared with the method according to drift-diffusion model in which impact ionization coefficient is calculated based only on the magnitude of the electric field, necessary time for calculation in such methods is significantly longer.